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1  Introduction 


In  imaging  systems  where  the  focal  plane  detector  array  is  not  sufficiently  dense,  so  as 
to  meet  the  Nyquist  criterion,  the  resulting  images  will  be  degraded  by  aliasing  effects. 
This  undersampling  is  a  common  problem  among  staring  infrared  imaging  systems.  This 
is  due  to  fabrication  complexities  and  quantum  efficiency  problems  associated  with  man¬ 
ufacturing  small  and  dense  infrared  focal  plane  arrays  (FPAs).  While  it  is  possible  to 
equip  such  systems  with  optics  to  properly  bandlimit  the  input,  this  generally  means 
employing  optics  with  a  very  small  instantaneous  field  of  view  (IFOV).  This  may  be 
highly  undesirable  in  some  applications.  Furthermore,  inexpensive  charge  coupled  device 
(CCD)  cameras  may  also  employ  detector  arrays  which  are  not  sufficiently  dense  for  the 
desired  optics.  Thus,  the  goal  of  this  work  is  to  obtain  high  resolution  images,  with 
reduced  aliasing,  from  such  systems.  In  addition,  we  wish  to  remove  the  blurring  effects 
of  the  finite  size  detectors  and  the  optics. 

One  way  to  overcome  the  undersampling  problem,  without  low  pass  filtering  the 
input  or  sacrificing  IFOV,  is  to  exploit  multiple  frames  from  an  image  sequence.  This  is 
possible  if  there  is  relative  motion  between  the  scene  and  the  FPA  during  image  sequence 
acquisition.  In  this  case,  a  unique  “look”  at  the  scene  (i.e.,  a  unique  set  of  samples)  may 
be  provided  by  each  frame.  The  desired  image  sequence  can  be  obtained  if  an  imager 
is  mounted  on  a  moving  platform,  such  as  an  aircraft.  We  refer  to  this  as  uncontrolled 
microscanning  [1].  It  is  also  possible  to  introduce  a  controlled  mirror  in  the  optical  path 
to  translate  the  scene  intensity  image  across  the  FPA.  This  is  referred  to  as  controlled 
microscanning  [1,2].  Here  we  focus  on  uncontrolled  microscanning  and  we  consider  both 
rotational  and  translational  motion  of  the  scene  relative  to  the  FPA.  The  key  to  the  high 
resolution  image  recovery  algorithm  is  accurate  knowledge  of  the  sub-pixel  translation 
and  rotation  of  each  frame.  If  these  parameters  are  unknown,  as  in  the  uncontrolled  case, 
they  must  be  estimated  from  the  observed  images.  Thus,  we  must  consider  both  image 
registration  and  high  resolution  image  reconstruction. 

Several  approaches  to  the  high  resolution  image  reconstruction  problem  have  been 
proposed  in  the  literature.  One  of  the  first  techniques  proposed  seeks  to  solve  for  an 
unaliased  discrete  spectrum  by  solving  a  set  of  frequency  domain  linear  equations  [3]. 
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This  approach  has  been  extended  to  treat  noise  in  [4]  and  blur  in  [5].  These  techniques, 
however,  do  not  address  rotation  which  may  be  present  in  many  applications.  Further¬ 
more,  we  have  found  that  they  tend  to  be  highly  sensitive  to  errors  in  the  registration 
parameters.  This  is  particularly  true  when  attempting  to  remove  blur  [5] . 

Another  approach  to  the  high  resolution  image  reconstruction  problem  uses  a  pro¬ 
jection  onto  convex  sets  (POCS)  algorithm  [6].  The  POCS  approach  has  been  extended 
to  treat  motion  blur  and  noise  in  [7,  8].  The  problem  has  also  been  approached  from  a 
statistical  estimation  framework.  Specifically,  a  maximum  a  posteriori  (MAP)  estimator 
is  developed  in  [9,  10]  which  is  an  extension  of  a  single  frame  image  expansion  algorithm 
proposed  in  [11].  A  statistical  estimation  approach  was  first  applied  to  forward  looking 
infrared  (FLIR)  data  in  [12].  In  particular,  a  maximum  likelihood  technique  using  the 
expectation  maximization  (EM)  algorithm  is  developed  which  seeks  to  jointly  estimate 
the  translational  shifts  and  a  high  resolution  image  [12].  In  addition,  a  joint  registration 
and  high  resolution  reconstruction  technique  using  MAP  estimation  is  presented  in  [13]. 
The  work  in  [12]  and  [13]  does  not  explicitly  treat  the  case  of  rotational  motion.  Another 
related  MAP  approach  can  be  found  in  [14]. 

The  reconstruction  algorithm  presented  here  can  be  viewed  as  an  extension  of  the 
basic  approach  presented  in  [15].  The  technique  in  [15]  seeks  to  minimize  a  specified 
cost  function  using  an  iterative  algorithm.  This  cost  function  is  the  total  squared  error 
between  the  observed  low  resolution  data  and  the  predicted  low  resolution  data.  The 
predicted  data  is  the  result  of  projecting  the  high  resolution  image  estimate  through  the 
observation  model.  Here  we  employ  the  same  registration  technique  used  in  [15]  and  we 
seek  to  minimize  a  related  cost  function.  However,  our  approach  includes  a  number  of 
important  and  fundamental  differences.  Some  of  the  novel  aspects  of  the  work  presented 
here  are  as  follows; 

•  The  observation  model  uses  accurate  information  about  the  optical  system  to  form 
a  realistic  point  spread  function  (PSF).  In  particular,  the  PSF  model  is  defined  to 
include  the  effects  of  the  optics  and  the  finite  size  detectors. 

•  The  cost  function  defined  here  includes  a  regularization  term.  This  extra  term  adds 
robustness,  particularly  when  only  a  small  number  of  frames  are  available  or  when 
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then  the  fidelity  of  the  data  is  low. 

•  The  unconstrained  minimization  of  the  cost  function  is  achieved  using  a  gradient 
descent  and  a  conjugate  gradient  technique.  The  conjugate  gradient  technique,  in 
particular,  provides  rapid  convergence. 

•  Finally,  we  study  the  application  of  high  resolution  image  reconstruction  to  a  real¬ 
time  infrared  imaging  system. 

The  organization  of  the  rest  of  the  paper  is  as  follows.  In  Section  2,  the  observa¬ 
tion  model  is  described.  Both  a  continuous  and  discrete  model  are  developed.  Section 
3  addresses  image  registration.  The  high  resolution  image  reconstruction  algorithm  is 
developed  in  Section  4.  In  particular,  the  regularized  cost  function  is  defined  and  two 
optimization  procedures  are  described.  Experimental  results  are  provided  in  Section  5. 
These  include  results  obtained  using  simulated  data  and  using  FLIR  images  acquired 
from  a  real-time  system.  Quantitative  error  analysis  is  provided  and  several  images  are 
shown  for  subjective  evaluation.  Finally,  some  conclusions  are  given  in  Section  6. 

2  Observation  Model 

In  this  section,  the  observation  model  is  presented.  This  model  is  the  basis  for  the  high 
resolution  reconstruction  algorithm  developed  in  Section  4.  We  begin  with  a  continuous 
model  which  closely  follows  the  physical  image  acquisition  process.  An  equivalent  discrete 
model  is  then  presented.  It  is  the  discrete  model  which  is  utilized  in  the  reconstruction 
algorithm.  We  conclude  this  section  with  the  characterization  of  the  system  point  spread 
function,  since  this  represents  a  key  element  in  the  observation  model. 

2.1  Continuous  Model 

A  block  diagram  of  the  continuous  observation  model  is  shown  in  Fig.  1.  The  true  scene 
intensity  image  is  defined  as  o(x,y).  The  motion  of  the  imager  that  occurs  between 
image  acquisitions  is  modeled  as  a  pure  rotation  and  translation  of  the  scene  intensity 
image.  For  a  moving  imager  and  a  stationary  scene  in  the  far  field,  this  is  a  fairly 


3 


Figure  1:  Continuous  observation  model  resembling  the  physical  process  of  image  acqui¬ 
sition. 

accurate  model,  since  occlusion  effects  and  perspective  changes  are  minimal.  Thus,  the 
A:’th  observed  frame  in  a  sequence  can  be  expressed  as 

Ok{x, y)  =  o{xcos6k  -  ysinOk  +  hk,  ycosBk  +  xsinOk  +  Vk),  (1) 

for  k  —  1, 2, 3, ...  ,p.  Note  that  6k  represents  the  rotation  of  the  A:’th  frame  about  the 
origin  (i.e.,  x  =  0,  y  =  0).  The  parameters  hk  and  Vk  represent  the  horizontal  and  vertical 
shift  associated  with  the  ^’th  frame. 

The  blurring  effect  of  the  optics  and  finite  detector  size  is  modeled  by  a  convolution 
operation  yielding 

dk{x,y)  =  Ok{x,y)*hc{x,y),  (2) 

where  hc{x,  y)  is  the  continuous  system  PSF.  More  will  be  said  about  the  PSF  in  Section 
2.3.  Finally,  the  blurred,  rotated  and  translated  image  is  sampled  below  the  Nyquist  rate 
and  corrupted  by  noise.  This  yields  the  k’th.  low  resolution  observed  frame 

yk{ni,n2)  =  Ok{niTi,n2T2)  +77fc(«i>^2),  (3) 

where  Ti  and  T2  are  the  horizontal  and  vertical  sample  spacings  and  %(ni,n2)  is  an 
additive  noise  term.  Let  the  dimensions  of  the  low  resolution  image  yfc(ni,  712)  be  N1XN2. 
These  data  in  lexicographical  notation  will  be  expressed  as  y*  =  [yk,iiyk,2T  •  ■  ,yk,M]'^, 
where  M  =  N1N2.  Finally,  let  the  full  set  of  p  observed  low  resolution  images  be  denoted 

y  =  [yi’>  y2 :  •  •  •  >  YpV  =  bi ,  yz,  •  •  • ,  ypMf-  (4) 

Thus,  all  observed  pixel  values  are  contained  within  y. 
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(Nyquist) 


Figure  2:  Equivalent  discrete  observation  model  illustrating  the  relationship  between  the 
ideally  sampled  image  z  and  the  observed  frames  y. 

2.2  Discrete  Model 

While  the  continuous  model  provides  insight  into  the  physical  process,  we  require  a 
discrete  observation  model  to  develop  the  high  resolution  reconstruction  algorithm.  That 
is,  we  need  a  model  relating  a  discrete  high  resolution  image  to  the  low  resolution  observed 
frames  y.  Figure  2  illustrates  such  a  discrete  model  which  is  equivalent  to  that  in  Fig. 
1.  The  difference  here  is  that  we  first  define  z{ni,n2)  to  be  an  intensity  image  sampled 
at  or  above  the  Nyquist  rate  with  no  blur  or  noise  degradation.  It  is  this  discrete  image 
we  wish  to  estimate  from  the  observed  frames.  Let  this  high  resolution  image  be  of  size 
LiNi  X  L2N2  =  N,  where  Li  and  L2  are  positive  integers.  More  will  be  said  about 
these  parameters  shortly.  In  later  analysis  it  will  be  convenient  to  express  this  image  in 
lexicographical  notation  as  the  vector  z  =  [zi,  2:2,  •  •  • »  ^nV- 

As  before,  this  image  is  rotated  by  6k  and  shifted  by  hk  and  Vk  producing  Zk(ni,  712). 
Here  we  will  define  the  shifts  hk  and  Vk  in  terms  of  low  resolution  pixel  spacings  for 
convenience.  Note  that  this  step  requires  interpolation  since  the  sampling  grid  changes 
in  the  geometric  transformation.  Theoretically,  one  could  use  ideal  interpolation  since 
z{ni,n2)  is  defined  to  be  sampled  above  the  Nyquist  rate.  However,  in  practice,  simpler 
interpolation  methods  such  as  nearest  neighbor  and  bilinear  interpolation  can  be  used 
[16].  For  large  values  of  Li  and  L2,  the  high  resolution  grid  is  so  dense  that  simple 
interpolation  methods  can  be  reasonably  accurate. 

Next,  the  system  PSF  is  accounted  for  yielding 

5(ni,n2)  =  Zk{nun2)  *  hd{nun2),  (5) 


5 


where  hd{ni,n2)  is  the  equivalent  discrete  system  PSF.  Note  that  the  blurring  is  per¬ 
formed  after  the  geometric  transformation.  If  the  PSF  is  circularly  symmetric,  then  the 
blurring  can  be  equivalently  introduced  prior  to  the  rotation.  This  saves  on  computations 
since  the  blurring  is  performed  only  once  in  this  case.  Finally,  the  transformed  image  is 
subsampled  down  to  the  resolution  of  the  observed  frames  yielding 

yk{ni,n2)  =  Zk{niLi,n2L2)  +  'nk{nun2).  (6) 

Aliasing  is  generally  introduced  in  this  final  step. 

This  discrete  model  can  be  rewritten  in  a  simple  form  where  the  low  resolution  pixels 
are  defined  as  a  weighted  sum  of  the  appropriate  high  resolution  pixels  with  additive  noise. 
This  generalized  form  can  account  for  the  PSF  blurring,  the  geometric  transformation  and 
the  subsampling  of  the  high  resolution  image.  Specifically,  the  observed  low  resolution 
pixels  from  frame  k  are  related  to  the  high  resolution  image  as  follows 

N 

yk,m  —  ^  ^  '^m,Tifik')  hk,  4-  T]k  (7) 

r=l 

for  771  =  1,2,...,M  and  k  =  1,2, ...,p.  The  weight  Wm,T{6k,hk,Vk)  represents  the 
“contribution”  of  the  r’th  high  resolution  pixel  to  the  m’th  low  resolution  observed  pixel 
of  the  fc’th  frame.  The  parameters  9k,  hk  and  Vk  represent  the  rotation,  horizontal  and 
vertical  translational  shifts,  respectively,  of  the  fc’th  frame  with  respect  to  some  reference 
on  the  high  resolution  grid.  The  term  r]k,m  in  (7)  represents  an  additive  noise  sample.  To 
further  compact  the  notation,  the  model  in  (7)  can  be  expressed  in  terms  of  the  entire 
set  of  low  resolution  pixels  as 

N 

ym  ~  ^  !  '^m,r^r  4"  (^) 

r=l 

for  m  =  1, 2, . . .  ,pM  and  where  Wm,r  is  simply  the  contribution  of  the  r’th  high  resolution 
pixel  in  z  to  the  m’th  low  resolution  pixel  in  y.  It  is  assumed  that  the  underlying 
scene,  z,  remains  constant  during  the  acquisition  of  the  multiple  low  resolution  frames. 
Furthermore,  we  assume  here  that  the  only  frame  to  frame  differences  in  the  weights 
result  from  rotation  and  translational  of  each  low  resolution  frame  relative  to  the  high 
resolution  grid. 
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(a)  (b) 


Figure  3:  Discrete  detector  model  showing  those  high  resolution  pixels  that  contribute 
to  a  low  resolution  pixel  for  two  different  registration  positions. 


A  simple  way  to  visualize  the  form  of  the  observation  model  in  (8)  is  to  consider 
only  the  blur  from  the  finite  detector  size.  This  scenario  is  illustrated  in  Fig.  3.  Here 
each  low  resolution  pixel  is  obtained  by  summing  the  high  resolution  pixels  within  the 
span  of  that  low  resolution  detector.  One  low  resolution  detector  in  Fig.  3(a)  is  shaded 
to  illustrate  this  point.  This  discrete  detector  model  simulates  the  integration  of  light 
intensity  that  falls  within  the  span  of  the  low  resolution  detector.  As  the  low  resolution 
grid  shifts  relative  to  the  fixed  high  resolution  grid,  as  in  Fig.  3(b),  a  different  set  of 
high  resolution  pixels  contribute  to  each  low  resolution  pixel.  This  yields  a  new  set  of 
linearly  independent  equations  from  (8).  Clearly,  some  type  of  interpolation  is  required 
for  any  non-integer  shift  on  the  high  resolution  grid  or  any  non-trivial  rotation.  This 
interpolation  can  be  accounted  for  by  modifying  the  weights  in  (8).  This  simple  detector 
model  can  give  good  results.  However,  a  more  realistic  PSF  model  is  described  in  the 
following  section. 
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2.3  System  Point  Spread  Function 

For  most  systems,  there  are  two  main  contributors  to  the  system  PSF.  The  primary 
contributor  is  generally  the  finite  detector  size  as  illustrated  in  Fig.  3.  This  effect  is 
spatially  invariant  for  a  uniform  detector  array.  The  second  contributor  is  the  optics. 
Here  we  assume  an  isoplanatic  model  for  the  optics  [17].  We  derive  and  use  a  theoretical 
PSF  because,  for  the  type  of  systems  considered  in  this  paper,  direct  measurement  of  an 
unaliased  system  PSF  is  not  possible. 

Let  us  begin  by  considering  a  system  with  a  uniform  detector  array.  The  effect  of  the 
integration  of  light  intensity  over  the  span  of  the  detectors  can  be  modeled  as  a  linear 
convolution  operation  with  a  PSF  determined  by  the  geometry  of  a  single  detector.  Let 
this  PSF  be  denoted  d{x,  y).  Applying  the  Fourier  transform  to  d{x,  y)  yields  the  effective 
continuous  frequency  response  resulting  from  the  detectors 

D{u,v)=J^T{d{x,y)},  (9) 

where  TT{-}  represents  the  continuous  Fourier  transform.  Next,  define  the  incoherent 
optical  transfer  function  (OTF)  of  the  optics  to  be  Hgiu,  v).  The  overall  system  OTF  is 
given  by  the  product  of  these,  yielding 

H{u,v)  =  D{u,v)Ho{u,v).  (10) 


The  overall  continuous  system  PSF  is  then  given  by 

h{x,y)=TT-^{H{u,v)}, 


(11) 


where  represents  the  inverse  Fourier  transform.  Finally,  the  impulse-invariant 

discrete  system  PSF  [18]  on  the  high  resolution  grid  is  obtained  by  sampling  the  contin¬ 
uous  PSF  such  that 


(12) 


L  /■„  „  \  ^  f  'n-iTi  n2T2\ 

hd(ni,n2)  -  j  . 

This  accurately  represents  the  continuous  blurring  when  the  effective  sampling  frequency 
Li/Ti  exceeds  two  times  the  horizontal  cutoff  frequency  of  H{u,v)  and  L2/T2  exceeds 
two  times  the  vertical  cutoff  frequency  [18]. 

Let  us  now  specifically  consider  a  system  with  uniform  rectangular  detectors.  An 
illustration  of  such  a  detector  array  with  critical  dimensions  labeled  is  provided  in  Fig. 
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Figure  4:  Uniform  detector  array  illustrating  critical  dimensions. 


4.  The  shaded  areas  represent  the  active  region  of  each  detector.  The  detector  model 
PSF  in  this  case  is  given  by 


y)  =  Irect  ^  ^ 

ab^^^Ka'b)  1  0  otherwise 


(13) 


Let  the  active  region  dimensions,  a  and  b,  be  measured  in  millimeters  (mm).  Thus,  the 
effective  continuous  frequency  response  resulting  from  the  detectors  is 


.  ,  /  ,  N  sin(7rau)sin(7r6?;) 

D(u,v)  =  smc{au,bv) - - , 


(14) 


where  u  and  v  are  the  horizontal  and  vertical  frequencies  measured  in  cycles/mm. 

The  incoherent  optical  transfer  function  (OTF)  of  diffraction-limited  optics  with  a 
circular  exit  pupil  can  be  found  [17]  as 


H,{u,v)  =  (  it)  -  t\l'-  -  (U'l  , 

[  0  otherwise 


(15) 


where  p  =  yjv?  +  v^.  The  parameters  pc  is  the  radial  system  cutoff  frequency  given  by 

1 


Pc 


A//#’ 


(16) 


where  //#  is  the  f-number  of  the  optics  and  A  is  the  wavelength  of  light  considered. 
Since  the  cutoff  of  Ho{u,  v)  is  pc,  so  is  the  cutoff  of  the  overall  system  H {u,  v).  Thus,  the 
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impulse-invariant  discrete  system  defined  in  (12)  will  accurately  model  the  continuous 
system  when  Lx  >  |'2pc^i]  and  L2  >  [2pcT2].  This  choice  of  Lx  and  L2  also  defines  a 
high  resolution  sampling  grid  at  or  above  the  Nyquist  rate  for  an  arbitrary  scene.  That 
is,  the  effective  high  resolution  sampling  rates  of  Lx/Tx  and  L2/T2  will  be  more  than 
twice  the  OTF  cutoff  frequency. 

Figure  5  shows  an  example  of  D{u,v),  Ho{u,v),  H{u,v)  and  hc{x,y)  for  a  particular 
imaging  system.  The  system  considered  is  the  FLIR  imager  used  to  collect  data  for  the 
experimental  results  presented  in  Section  5.  The  FLIR  camera  uses  a  128  x  128  Amber 
AE-4128  infrared  FPA.  The  FPA  is  composed  of  Indium-Antimonide  (InSb)  detectors 
with  a  response  in  the  3/xm  -  5//m  wavelength  band.  This  system  has  square  detectors 
of  size  a  =  b  =  .040mm.  The  imager  is  equipped  with  100mm  //3  optics.  The  center 
wavelength,  A  =  .004mm,  is  used  in  the  OTF  calculation.  Figure  5(a)  shows  the  effective 
modulation  transfer  function  (MTF)  of  the  detectors,  |Z1(m,  u)|.  The  diffraction-limited 
OTF  for  the  optics,  Ho{u,v),  is  shown  in  Fig.  5(b).  Note  that  the  cutoff  frequency  is 
83.3  cycles/mm.  The  overall  system  MTF,  \H{u,v)\,  is  plotted  in  Fig.  5(c).  Finally,  the 
continuous  system  PSF,  hc{x,y),  is  plotted  in  Fig.  5(d). 

The  detector  spacing  on  the  Amber  FPA  is  Tx=T2  =  .050mm,  yielding  a  sampling 
frequency  of  20  cycles/mm  in  both  directions.  Thus,  the  effective  sampling  rate  must  be 
increased  by  a  factor  of  8.33  to  eliminate  aliasing  entirely  for  an  arbitrary  scene.  This 
would  require  that  we  select  Lx,  L2>  9.  In  practice,  we  find  that  good  results  can  be 
obtained  with  smaller  values  of  Lx  and  L2- 

3  Image  Registration 

In  most  applications,  the  registration  parameters  in  the  observation  model,  9^,  hk  and 
Vk,  will  not  be  known  a  priori.  Thus,  they  must  be  estimated  from  the  observed  image 
sequence.  Accurate  sub-pixel  registration  is  the  key  to  the  success  of  the  high  resolution 
image  reconstruction  algorithm.  A  number  of  image  registration  techniques  have  been 
proposed  in  the  literature.  We  have  found  that  a  practical  and  effective  method  of 
estimating  the  sub-pixel  translation  and  rotation  is  using  the  technique  in  [15],  which  is 
based  on  that  in  [19].  For  convenience,  this  algorithm  is  presented  here  using  the  current 
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Continuous  Detector  MTF 


Continuous  OTF  of  Optics 


(a) 


(b) 


Continuous  Overall  System  MTF 


Continuous  System  PSF 


Figure  5:  (a)  Effective  MTF  of  the  detectors  in  the  FLIR  imager  (b)  diffraction-limited 
OTF  of  the  optics  (c)  overall  system  MTF  (d)  overall  continuous  system  PSF. 
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notation. 


To  begin,  define  the  first  observed  frame  to  be  the  reference  frame  and  without  loss 
of  generality  let  0i  =  hi  =  vi  =  0.  According  to  our  model, 

Ok {x,  y)  =  Oi  {xcosOk  -  ysin^fe  +  hk,  ycosBk  +  rrsin^fc  +  ,  (17) 

for  k  =  2, 3, . . .  ,p.  Note  that  this  assumes  that  the  center  of  rotation  is  at  the  origin 
(i.e.,  X  =  0,  y  =  0).  This  is  not  restrictive  however,  since  we  allow  any  shift  hk  and  Vk- 
If  the  PSF  blur  is  approximately  circularly  symmetric,  then 

Ok  (x,  y)  «  di  {xcosOk  -  ysvcidk  +  hk,  ycosdk  +  xsinOk  +  Vk).  (18) 

For  very  small  values  of  6k,  we  can  make  the  following  approximations:  sin^*  ^  9k  and 
cosOk  ~  1.  Using  these  yields 

Ok{x, y)  ~  oi{^  -  y^k  -khk,y  +  xOk  +  Vk).  (19) 

Now  we  use  the  first  three  terms  of  the  Taylor  series  expansion  as  an  approximation  for 
the  right  side  in  (19).  This  yields 

dk{x,y)  «  di{x,y)  +  {hk  -  y9k)gx{^,y)  +  {vk  +  x6k)gy{x,y),  (20) 

where  gx{x,y)  =  and  gy{x,y)  = 

In  light  of  the  relationship  expressed  in  (20),  we  define  the  least  squares  estimates  for 
the  registration  parameters  as  follows 

arg  min 

6kj  hk,  Vk  —  6k,  hk,  Vk  Ek{6k,  hk,  Vk),  (^1) 


where 

Ek{6k,  hk,  Vk)=  {dk{x,y) -di{x,y)  -  {hk-y6k)gx{x,y)  -  {vk  +  x6k)gy{x,y))'^ . 

{x,y)es 

(22) 

Note  that  S  represents  the  grid  of  points  in  the  space,  defined  by  x  and  y,  at  which 
we  have  discrete  samples.  Rewriting  this  error  in  terms  of  the  observed  images  yields 

Ek{6k,  hk,  Vk)  =  Y,  -  2/1  (’^)  “  (^^  ~  n2T26k)gx{n)  -  {vk  +  niTi6k)gy{n))‘^ , 

nQAf 

(23) 
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where  n  =  [ni,  712]  and  M  is  the  set  of  indices  on  the  low  resolution  discrete  grid  for  which 
we  have  observations.  Note  that  the  center  of  rotation  on  the  discrete  grid  is  assumed  to 
be  at  m  =  712  =  0.  The  functions  ^^(n)  and  pj,(n)  are  discrete  estimates  of  gx{x,y)  and 
gy{x,y),  respectively,  at  location  x  =  riiTi  and  y  =  n2T2.  These  can  be  computed  using 
scaled  Prewitt  operators  [16],  for  example. 

To  solve  the  minimization  problem  in  (21),  we  begin  by  differentiating  E{9k,  hk,  u*) 
with  respect  to  6k,  hk  and  Vk  and  set  the  derivatives  equal  to  zero.  This  yields  the 
following  three  equations 


{hkglin)  +  Vkgx{n)gy{n)  +  9kg{n)ga:{n))  =  yfe(n)5x(n) 


n£Af 


and 


where 


and 


(Mi(n)^j,(n)  +  UfcpJ(n)  +  9kg{n)gy{n)'j  =  ^  ^^(11)52,(11), 
neAf  nGAT 

{hkg{n)g^{n)  +  Vkg{n)gy{n)  +  9kf{n)')  =  yk{n)g{n), 
neAf  neAf 


g{n)  =  niTigy{n)  -  n2T2^x(n) 


Vkin)  =  yk{n)  -  yi{n). 

We  then  simultaneously  solve  these  expressions.  To  do  so  let 


MRk  =  Vk, 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


where 


M  = 


EneAr5x(n)  9xin)gy{n)  Zne^^ 9{^)9x{n) 

^ne^f  9xi^)9y{^)  SneAC  ^(^)5j/(^) 


? 


Rk  —  \hk,'^k,9k\  ,  and 


14  = 


EneAryfc(n)5x(n)  ' 
EneJV  yfc(n)P2/(n) 
EneAr2/fc(n)5(n)  _ 


Finally,  the  estimated  registration  vector,  Rk  =  [hk,  Vk,  9k\^,  can  be  computed  as 


(30) 


(31) 


Rk  =  M-^Vk. 


(32) 
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Table  1:  Iterative  registration  algorithm  for  yfc(n),  where  k  =  2,3, . . .  ,p. 


step  1:  Compute  5j;(n),  ^y(n)  and  p(n)  from  t/i(n)  and  form  M  matrix. 

step  2:  Let  n  =  0,  y°(n)  =  and  =  [0, 0, 0]^. 

step  3:  Compute  T4  using  yfc(n)  =  yfc(n)  —  yi(n)  and  then  let  n  =  n  +  1. 

step  4;  Compute  =  M~^Vk  +  Rk~^- 

step  5:  if  \\R^  -  .Rfc~^||/||^fc"^||  <  T,  let  Rk  =  Rl  and  stop. 

step  6:  Resample  yfc(n)  towards  yi(n)  according  to  R^  to  create  yjt(n)  and 
go  to  step  3. 


To  obtain  shifts  in  terms  of  low  resolution  pixel  spacings  (rather  than  in  mm),  we  set 
Ti  =  T2  =  1  in  (27). 

Because  of  the  assumptions  made,  this  technique  is  only  accurate  for  small  shifts  and 
rotations.  To  treat  the  case  where  larger  values  are  expected,  we  follow  the  iterative 
method  used  in  [15]  and  [19].  To  do  so,  the  initial  registration  parameters  are  estimated 
according  to  (32).  Next  yfc(n)  is  shifted  and  rotated  according  to  the  registration  param¬ 
eter  estimates  so  as  to  more  closely  match  yi(n).  This  modified  image  is  then  registered 
to  yi(n).  The  process  continues  whereby  yA;(ii)  is  continually  modified  until  the  registra¬ 
tion  estimates  become  suflBciently  small.  The  final  registration  estimate  is  obtained  by 
summing  all  of  the  “partial”  estimates.  The  iterative  registration  procedure  for  yfc(n), 
where  A:  =  2, 3, . . .  ,p,  is  summarized  in  Table  1. 

Because  the  three  parameters  are  well  overdetermined  by  the  data,  this  least  squares 
estimate  is  generally  accurate.  We  find  that  the  main  source  of  error  lies  in  the  resam¬ 
pling  of  yfc(n)  (step  6  in  Table  1),  since  this  requires  interpolation  on  the  low  resolution 
aliased  grid.  Some  error  is  also  introduced  in  the  discrete  gradient  estimates.  In  general, 
however,  the  algorithm  appears  to  provide  sufficiently  accurate  results  for  our  current 
application.  In  cases  where  the  aliasing  is  more  severe,  a  joint  registration  and  recon- 
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struction  algorithm,  like  that  described  in  [12],  [14],  and  [13],  may  be  advantageous. 


4  High  Resolution  Image  Reconstruction 


With  estimates  of  the  registration  parameters,  the  observation  model  can  be  completely 
specified.  In  light  of  the  observation  model  in  (8),  we  define  the  high  resolution  image 
estimate  to  be 


where 


and  for  m  =  1, 2, . . .  ,pM  are  the  observed  pixel  values.  Clearly,  the  cost  function  in 
(34)  balances  two  types  of  errors.  The  first  term  on  the  right  hand  side  is  minimized  when 
z,  projected  through  the  observation  model,  matches  the  observed  data.  However,  direct 
minimization  of  this  term  can  lead  to  excessive  noise  magnification  due  to  the  ill-posed 
nature  of  the  inverse  problem.  Thus,  the  second  term  serves  as  a  regularization  operator. 
The  parameters  ai^j  are  generally  selected  so  that  the  regularization  term  is  minimized 
when  z  is  smooth.  Here  we  select  aij  to  be  1/4  only  for  those  four  values  of  j  which 
correspond  to  immediate  spatial  neighbors  of  Zi.  The  others  are  set  to  zero.  The  weight 
of  these  competing  “forces”  in  the  cost  function  is  controlled  by  the  “tuning”  parameter 
A.  Larger  values  of  A  will  generally  lead  to  a  smoother  solution.  This  is  useful  when  only 
a  small  number  of  frames  is  available  or  the  fidelity  of  the  observed  data  is  low.  It  is  also 
possible  to  make  A  spatially  adaptive  in  a  fashion  similar  to  that  in  [20,  21].  Finally,  note 
that  the  estimate  defined  in  (33)  and  (34)  can  be  viewed  as  MAP  estimate  in  the  case  of 
Gaussian  noise  and  where  z  is  a  realization  of  a  Gauss-Markov  random  process  [13]. 

Next  we  consider  two  unconstrained  optimization  techniques  for  minimizing  the  cost 
function  in  (34).  First  we  consider  a  gradient  descent  approach  and  then  we  present  a 
conjugate  gradient  method. 
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4.1  Gradient  Descent  Optimization 


To  derive  the  gradient  descent  update  procedure  for  the  image  estimate,  we  begin  by  dif¬ 
ferentiating  (34)  with  respect  to  some  pixel  Zk  ioi  k  =  1,2, ...  ,N.  This  partial  derivative 
is  given  by 


9k{z)  = 


dC{z) 

dzk 


pM 

m=l 


r  N 


Vm  | 

Lr— 1 

N  N 

+  A  j  2:*  —  ^  akjZj  -t-  ^  ai^k 

j=l  2=1 


(35) 


N 

E 

J  =  1 


OiijZj  Zi 


The  iterative  procedure  begins  with  an  initial  estimate  of  the  high  resolution  image  z^. 
A  relatively  simple  starting  point  can  be  obtained  by  interpolating  the  first  frame  using 
bilinear  or  bicubic  interpolation.  The  gradient  descent  update  for  each  pixel  estimate  is 


-71+1 


(36) 


for  n  =  0, 1, 2, . . .  and  k  =  l,2,...,N.  Alternatively,  the  update  can  be  written  as 


^n+l  ^ 


(37) 


where 

‘  i?i(z")  " 

n  ^2(z”) 

s  =  : 

.  5iv(z”)  . 

The  parameter  in  (36)  and  (37)  represents  the  step  size  at  the  n’th  iteration.  This 
parameter  must  be  selected  to  be  small  enough  to  prevent  divergence  and  large  enough 
to  provide  convergence  in  a  reasonable  number  of  iterations.  The  optimal  step  size  can 
be  calculated  by  minimizing 

(7(2"+^)  =  (7(2"  -  eV)  (39) 


with  respect  to  £■".  To  do  so  we  begin  by  writing  using  (34)  and  (36).  Next  we 

differentiate  this  with  respect  to  e"  and  set  the  derivative  equal  to  zero.  Solving  for  e" 
yields,  after  some  manipulation, 

„  7..  (£^,1  -  s/„)  +  A  a  (if  -  £f.i  a,jz]) 


a"i75.  +  AE£iffi 


N  7;2 

i 


(40) 
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Table  2;  Proposed  gradient  descent  iterative  estimation  algorithm. 


step  1:  Begin  at  n  =  0  with  initial  estimate  z°  being  the  interpolated  low 
resolution  frame  1. 

step  2:  Compute  the  gradient  ^fc(z")  given  in  (35)  for  k  =  1,2, . . .  ,N,  yielding 

gn_ 

step  3:  Compute  the  optimal  step  size  e”  using  (40). 
step  4:  Let  z"'*’^  =  z"  - 

step  5:  If  ||z"+’-  -  z”||/||z"||  <  T,  let  z  =  z"+^  and  stop, 
step  6:  Let  n  =  n  +  1  and  go  to  step  2. 


where  ^ 

7m  =  l^^m,r5r(z")  (41) 

r=l 

is  the  gradient  projected  through  the  model,  and 

gi  =  Qi (z”)  -  5]  otijQj {tT)  (42) 

is  c?i(z”)  minus  the  weighted  sum  of  its  “neighbors.”  This  iteration  continues  until  the 
cost  function  stabilizes  or  |lz"+^  -  z”l|/|lz"||  <  T,  where  T  is  a  specified  threshold  value. 
A  summary  of  the  gradient  descent  optimization  procedure  is  provided  in  Table  2. 

4.2  Conjugate  Gradient  Optimization 

In  this  section,  we  describe  a  conjugate  gradient  optimization  procedure  for  minimizing 
the  cost  function  in  (34).  In  particular,  we  employ  the  Fletcher-Reeves  method  [22].  We 
later  show  that  with  little  additional  computational  complexity,  faster  convergence  can 
be  achieved  using  this  method  compared  to  gradient  descent. 

The  basic  conjugate  gradient  image  update  is  given  by 

z^1  =  4"  +  £"4(z"),  (43) 
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for  n  =  0,1,2,...  and  k  =  1,2,  Here  dfc(z")  is  the  conjugate  gradient  term. 

Alternatively,  the  update  can  be  written  as 

z"+i  =  z"H-eM",  (44) 


where 


d" 


di(z”)  ' 

d2(z") 


L  dN{z^)  J 


(45) 


As  before,  the  parameter  e”  is  the  step  size  at  the  n’th  iteration.  The  optimal  step  size 
can  be  calculated  by  minimizing  C'(z"‘''^)  =  (^(z"  +  e"d”)  with  respect  to  £■”.  This  yields 


6^^  = 


<^m  {Er=l  -  t/m)  +  A  dj  (if  -  Ef=l  C^i,jZj) 


(46) 


where 


N 

r=l 


and 


N 


di  =  dj(z")  -'^aijdj{z^). 
j=i 


The  conjugate  gradient  vector  is  initially  set  to  be 


(47) 

(48) 


d°  = 


(49) 


This  is  updated  according  to 


jjn+l  ^  _gn+l  ^ 


_L 


(50) 


where 


13^  = 


(g7i+l)Tgn+l 


(51) 


(gn)Tgn 

Again,  the  iterations  continue  until  the  estimate  converges.  A  summary  of  the  conjugate 
gradient  optimization  procedure  is  given  in  Table  3. 
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Table  3:  Proposed  conjugate  gradient  iterative  optimization  algorithm. 


step  1:  Begin  at  n  =  0  with  initial  estimate  z°  being  the  interpolated  low 
resolution  frame  1. 

step  2:  Compute  g°  and  initialize  conjugate  gradient  vector  as  d°  =  —  g°. 
step  3:  Compute  the  optimal  step  size  e"  using  (46). 
step  4:  Let  =  z"  + 

step  5:  If  llz”"''^  —  z”||/||z”||  <  T,  let  z  =  and  stop. 

step  6:  Compute  g""^^  and  let  d”"*"^  =  — g""^^  +  /?”d”,  where  /?”  =  ^^(gnjrgn'*'^  • 

step  7:  Let  n  =  n  +  1  and  go  to  step  3. 


5  Experimental  Results 

In  this  section  a  number  of  experimental  results  are  presented  in  order  to  demonstrate  the 
performance  of  the  proposed  algorithm.  The  first  set  of  experiments  involve  simulated 
data  derived  from  a  single  broad  band  visible  image.  The  second  set  of  results  use  data 
obtained  from  a  real-time  FLIR  imaging  system. 

5.1  Simulated  Imagery 

Here  we  use  simulated  data  in  order  to  evaluate  the  algorithms  quantitatively.  In  par¬ 
ticular,  16  low  resolution  images  are  generated  by  rotating,  translating,  blurring  and 
subsampling  an  “ideal”  image.  The  rotation  and  translation  parameters  for  each  frame 
have  been  selected  arbitrarily  and  are  shown  in  Fig.  6.  The  original  image  is  of  size 
250  X  250  and  the  down  sampling  factors  are  Li  =  L2  =  5.  The  blurring  function  is 
a  5  X  5  moving  average  filter,  which  simulates  the  low  resolution  detector  effects.  Fi¬ 
nally,  additive  Gaussian  noise  of  variance  =  10  is  introduced  in  each  frame.  The 
estimated  registration  parameters,  using  the  algorithm  in  Table  i,  are  also  shown  in  Fig. 
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Registraion  Parameters 


Figure  6:  Actual  and  estimated  registration  parameters  for  each  of  the  16  simulated  low 
resolution  frames.  The  shifts  hk  and  Vk  are  measured  in  terms  of  low  resolution  pixel 
spacings. 

6.  The  mean  absolute  error  (MAE)  between  the  actual  and  estimated  translational  shifts 
is  0.0414  low  resolution  pixel  spacings  and  the  MAE  for  the  rotation  parameters  is  0.0250 
degrees. 

The  original  8  bit  grayscale  image  “Aerial”  is  shown  in  Fig.  7(a).  One  typical  frame 
of  the  simulated  low-resolution  noisy  data  is  shown  in  Fig.  7 (b) .  The  multi- frame  image 
estimate,  using  the  estimated  registration  parameters,  is  shown  in  Fig.  7(c).  The  initial 
image  estimate  used  here  is  a  bilinearly  interpolated  version  of  the  first  frame.  Twenty 
iterations  of  the  conjugate  gradient  optimization  have  been  performed  with  A  =  0.1.  We 
find  that  the  algorithm  is  not  highly  sensitive  to  the  choice  of  A.  For  comparison  with  the 
multi-frame  estimate,  a  bicubic  interpolation  of  the  single  frame  in  Fig.  7(b)  is  shown  in 
Fig.  7(d). 

To  show  the  convergence  behavior  of  the  gradient  descent  and  conjugate  gradient 
procedures,  the  values  of  the  cost  function  in  (34)  are  plotted  in  Fig.  8  versus  iteration 
number  for  each.  Note  that  the  conjugate  gradient  method  shows  superior  convergence 
speed.  The  MAE  between  the  multi-frame  estimate  and  the  “ideal”  image  is  plotted 
in  Fig.  9  as  a  function  of  the  number  of  frames  used.  For  comparison,  the  MAEs 
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Original  image  Aerial 


50  100  150  200  250 


(a) 


Noisy  Low-Resolution  Observed  Frame  1 


10  20  30  40  50 


(b) 

Figure  7:  (a)  Original  image  “Aerial”  (b)  simulated  low  resolution  frame  1  where  = 
L2  =  5  and  =  10  (c)  multi-frame  image  estimate  using  16  frames  with  A  =  0.1  (d) 
bicubic  interpolation  of  frame  1. 
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Cost  Function  vs.  Iterations 


Figure  8:  Convergence  behavior  of  the  gradient  descent  and  conjugate  gradient  proce¬ 
dures  using  the  16  frames  of  simulated  data.  The  noise  variance  is  <t^  =  10  and  A  =  0.1. 

of  the  first  frame  bilinearly  and  bicubically  interpolated  are  also  shown.  With  only 
one  frame,  the  proposed  algorithm’s  performance  is  only  slightly  better  than  that  of  the 
bicubic  interpolator.  However,  with  additional  frames,  the  estimate  becomes  significantly 
improved  with  respect  to  the  single  frame  interpolators.  Similar  results  are  observed  with 
mean  squared  error. 

5.2  FLIR  Imagery 

Now  we  consider  applying  the  multi-frame  algorithm  to  FLIR  data  from  the  imager 
described  in  Section  2.3.  We  have  chosen  to  perform  a  reconstruction  with  L1  —  L2  —  5, 
although  this  resolution  may  not  be  sufficient  to  avoid  aliasing  effects  entirely.  The 
theoretical  discrete  PSF  of  the  FLIR  system  on  the  high  resolution  grid,  given  by  (12), 
is  shown  in  Fig.  10. 

The  multi-frame  algorithm  is  tested  using  20  frames  of  the  FLIR  imagery.  The  frames 
have  been  acquired  at  a  60  frame  per  second  rate.  The  rotation  and  translation  is  in¬ 
troduced  by  arbitrarily  manipulating  the  imager  during  acquisition.  One  typical  original 
resolution  frame  is  shown  in  Fig.  11(a).  This  80  x  80  image  represents  a  region  of  interest 
from  the  full  128  x  128  Amber  array.  The  scene  contains  a  number  of  small  power  boats 
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MAE  vs.  Number  of  Frames  Used 


Figure  9:  Mean  absolute  error  for  the  multi-frame  estimator  using  different  numbers  of 
frames.  The  noise  variance  is  <7^  =  10  and  A  =  0.1. 

and  trailers  on  a  gravel  parking  lot  with  a  fence  in  the  foreground.  The  multi-frame 
estimate  is  shown  in  Fig.  11(b)  for  A  =  0.1.  Again,  the  initial  image  estimate  is  obtained 
by  bilinearly  interpolating  the  first  frame  and  10  iterations  of  the  conjugate  gradient 
procedure  have  been  performed.  For  comparison,  a  bilinearly  interpolated  single  frame 
is  shown  in  Fig.  11(c).  A  bicubically  interpolated  version  of  the  single  frame  is  shown 
in  Fig.  11(d).  The  multi-frame  reconstruction  appears  to  show  significantly  improved 
image  detail.  In  addition,  note  that  the  aliasing  artifacts  on  the  diagonal  beam  of  the 
gate  in  the  foreground  of  Figs.  11(a),  11(c)  and  11(d)  are  virtually  eliminated  in  the 
multi-frame  estimate.  The  estimated  registration  parameters  for  the  20  observed  frames 
are  shown  in  Fig.  12. 

Finally,  to  illustrate  the  convergence  behavior  of  the  algorithms  using  the  FLIR  data, 
the  cost  function  is  plotted  in  Fig.  13  versus  iteration  number  for  both  the  gradient 
descent  and  the  conjugate  gradient  optimization  methods.  Again,  the  conjugate  gradient 
algorithm  exhibits  slightly  faster  convergence. 
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Discrete  System  PSF 


Figure  10:  Theoretical  discrete  system  PSF  of  the  FLIR  imager  for  =  L2  =  5. 

6  Conclusions 

Aliasing  reduction  and  resolution  enhancement  can  be  achieved  by  exploiting  multiple 
frames  which  are  rotated  and/or  translated  with  respect  to  one  another.  This  is  possible 
because  each  frame  offers  a  unique  set  of  discrete  samples.  For  an  imager  mounted 
on  a  moving  platform,  such  as  an  aircraft,  the  desired  image  sequence  may  arise  from 
natural  line-of-sight  jitter  and  rotation  of  the  platform.  With  this  in  mind,  it  may  then 
be  possible  to  relax  image  stabilization  requirements  in  some  applications  and  obtain 
improved  resolution  images  through  the  proposed  algorithm.  The  tradeoff  here  is  that 
of  temporal  resolution  for  spatial  resolution.  That  is,  multiple  low  resolution  frames  are 
required  to  form  each  new  high  resolution  image. 

The  key  to  the  success  of  the  algorithm  is  having  an  accurate  observation  model.  This 
includes  the  image  registration  parameters  and  the  system  PSF.  The  observation  model 
proposed  here  includes  accurate  information  about  the  optical  system.  We  believe  that 
this  provides  a  reasonably  accurate  system  PSF.  A  regularized  cost  function  defines  the 
image  estimate.  Minimization  of  the  cost  function  is  performed  using  either  a  gradient 
descent  or  conjugate  gradient  technique. 

The  quantitative  results  obtained  show  that  the  multi-frame  image  estimates  have  sig¬ 
nificantly  lower  error  than  estimates  formed  by  single  frame  interpolation.  Furthermore, 
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Registralon  Parameters 


h_k 


Figure  12:  Estimated  registration  parameters  for  the  20  frames  acquired  with  the  FLIR 
imager. 

we  believe  that  the  FLIR  results  show  that  the  multi-frame  estimate  has  significantly- 
improved  image  detail.  In  particular,  edges  and  fine  structure  emerge  in  the  multi-frame 
reconstruction  that  are  not  visible  in  the  low  resolution  data.  Because  these  features 
offer  important  visual  cues,  we  believe  that  the  utility  of  the  processed  image  is  greatly 
enhanced. 
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Cost  Function  vs.  Iterations 
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